Wind driven 3D Navier-Stokes circulation in the Atlantic 
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Abstract 

A finite element method for the numerical solution of the anisotropic Navier-Stokes equations 
in shallow domain is presented. This method take into account aspect ratio in the hydrostatic 
approximation of the Navier-Stokes equations [H [51 A projection method [71 [TJ] is used 
for the time discretization. The linear systems are solved via a some preconditioned conjugate 
algorithm, well adapted to massively parallel computers [17 1 116 1 [T8 ] . Some results are presented 
for the wind driven water circulation in the North Atlantic. 

Keywords Navier-Stokes equations, shallow water, projection methods, parallel numerical linear 
algebra, preconditioned conjugate gradient. 

1 Introduction 

The numerical simulation of wind-driven currents in the ocean is a common and current subject. 
Most of the used models are based on the hydrostatic approximation and anisotropic viscosity. A 
Galerkin method for the primitive equations in shallow domains is presented in [ 101 tH] • Prismatic 
finite elements are used, and numerical example are presented for the English Channel. In [6j, a 
baroclinic circulation model for the North Atlantic is presented. Spherical-polar coordinated are 
used . A study of the interaction between a baroclinic and a barotropic model is given in |12| . It 
is based on the planetary geostrophic equations. Finally in [20] a barotropic model is used for the 
study of the wind-driven circulation in an elongated basin. Most of the anisotropic models used in 
oceanography can be mathematically justified; see e.g. [H |2l [T9l [3]. 

Since the publication of the M.S. Lozier paper: Deconstructing the Conveyor Belt [9], it is of 
importance to have a good description of the wind-driven circulation in oceans. This paper is 
devoted to the numerical simulation of such a circulation in the North Atlantic, with a particular 
attention to the 3D aspects; including down and up-welling. 

In this paper, the full system of Navier-Stokes equations is considered; it is organized as follows. 
In section [2| we recall the setting of the anisotropic Navier-Stokes equations in the context of thin 
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domains. Then a description of a projection method is presented in section [3| A weak formulation 
and a finite element discretization are given in sections |4] and [5} We end this paper with section [6] 
where some numerical simulations obtained with our method are presented. Some wind driven 3D 
Navier-Stokes circulation in the North Atlantic is reported. 

2 Anisotropic Navier-Stokes equations 

Water flows in oceanography and limnology are governed by the Navier-Stokes equations. For some 
numerical simulations, asymptotic models are in current use (see |13t I21j). These are all based on 
the fact that the horizontal dimensions of the considered domain are much larger than the vertical 
one. Let d be the horizontal width, and h be the depth. The simplest model using the fact that 

h 

is very small is the hydrostatic model. In this model, we take care of turbulence effects by setting 
an anisotropic viscosity, much smaller in the vertical direction than in the horizontal one (see ^13j). 
Let ri = C be the domain defined by 

0. = [x = (xi,X2,X3) e M^, (xi,X2) G Ts, -h{xi,X2) < X3 < 0} 

where Tg is the surface of the domain and /i : — )■ M is its depth. The bottom of the domain is 
defined by Tb = dn\ F, (Fig. 

d 



Wind 




Figure 1 : Illustration of the shallow domain 

It is assumed that the water motion is generated by horizontal tension, induced by some wind on 
the surface F^. This motion is driven by the following anisotropic Navier-Stokes equations and is 
influenced by the Coriolis force: 

ut + {u \ V) u = AyU - 2uj Au -Vp in t > 0, (1) 
div n = inn, t> 0, (2) 
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with the following boundary and initial conditions. The bottom is subdivided into two parts Tq 
with strictly positive measure, and Fi = \ Fq, and 



u = on Fo, t > 0, (3) 

(u I n) = 0, (fj -n I r) = onFi, t>0, (4) 

where a = aij is the stress tensor, with a = —pi + Vti ■ v + v ■ V«*, n is the unit outward normal 
vector to the boundary, and r is any tangential vector. Note that the part Fi of the bottom 
corresponds to artificial limits of the domain. 
On the surface F^, tension conditions are considered 

'^^dx3 ^ '^^dxs ^ "3 = on F^, t > 0. (5) 
Finally the initial condition is 

u{-,t = 0) = inn. (6) 

The following notations have been used: u = (^1,^2,^3) is the fluid velocity, a; = {0,0, uj^) is the 
angular velocity of the Earth (projected onto the vertical in local coordinates), p is the pressure, 
u = diag(i/i, f2, fs) is the turbulent viscosity diagonal tensor, 9i and 62 are the tensions induced 
by the wind and 

j=l ^ 

Let us do the following change of variables and functions 

xi := xi, X2--X2, xs'-xs/e 
vi := ui, V2 ■= U2, v^:=u^/e 

With this scale change we get 

+ Vvi) - l^^evi -fv2 + g^^ = in t > 0, 
+ (t; I Vv2) - A,.V2 + /^i + ^ = mn,t>0, 

e''{^ + {v\Vv3)-A,eVs} + ^^ = mn,t>0, 

div V = inU, t> 0, 

vs- — = €01, 1^37— = €6*2, ^3 = on Fs, t>0, 
0x3 0x3 

V = on Fo, i > 0, 

{v \n) = 0, (o-f • n I r) = on Fi, t > 0, 

with = diag(i^i, U2, ^'3/e^), and / = 20J3. 
Set 

J^l = Al, f2 = A2, 1^3 = ^ A3, 
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Q, = e-^ei i = l,2, 
A = (Ai, A2, A3). 

When e — >• 0, this problem becomes the hydrostatic approximation of Navier-Stokes equations. 

^ + iv\ ^vi) - ^Wi -fv2 + ^^ = in t > 0, (7) 
-^ + iv\ V^^2) - AxV2 + fvi + Q^^ = in i > 0, (8) 

^ = in t > 0, (9) 

9X3 

diYV = in 17, t > 0. (10) 

with 

^1 = ""2 = "^3 • "-3 = on To, t > 0, (11) 
(VH Inn) = 0, {an -nn \th) = onTi, t>0, (12) 
where wh denotes the horizontal components of w, 

h^ = Qi, A37^ = e2, ^^3 = onr„t>o, (13) 
0x3 0x3 

vi{-^t = Q) = V2{-,t = Q) = Q in J]. (14) 

This development shows that it is natural to use an anisotropic viscosity diagonal tensor with the 
third component of the order of e^, compared with the others. 

3 Projection method 

In order to solve equations Q-Q, a velocity-correction projection method [71 [T51 E] is used. Let 
us recall this method in our case. Set vP = u{0),p^ = p{0), and let u^, and p^ be approximations 
of u{6t) and p{5t). A BDF2 scheme is used for the time discretization. For k > 1, we look for 



, and u^~^^ such that the prediction n'^^^ of the velocity is a solution of 
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5t\2 



^ - A^u''-^^ = -(u''\ v) u'' -2loAu'' - Vp^ in Q, (15) 

= on To, (16) 
n'=+^ \n) = on Ti, (17) 



= onr, (18) 

8X3 

= 0',^\ onr, (19) 

0x3 

^3+^ = on r^, (20) 
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then the velocity u'^"'"^, and the pressure satisfy 



26t 



,k+l 





= 


in 17, 




= 


in 17, 


n\ 


= 


on do. 



If q^~^^ = — ^ from equations (21 ) and (22 ) we have 

3 



dn 



26t 
on dn 



divtt^+^ in Q, 



(21) 
(22) 
(23) 

(24) 
(25) 



qk+l ^ pk^ ^k+l ^ ~k+l ^ St-Vq''^^. 
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If q^~^^ is known, then p^+i 

Note that this method consists in computing a non divergence-free prediction u^~^^ of the velocity, 
and then obtain u^~^^ as a divergence-free vector field via a pressure correction. 



4 Weak formulation 



A week formulation of problems O5P-(^20p and (24)-(25) is the following. Define 



V = {ip£H^{n); (^ = Oonro}, 
W = G (u I n) = Oonr, UTi} , 
U = {v £ H^{Qf; {v\n) = Oondn} . 

For 6*1, 02 G H~^/'^{rs), if /c > 0, assume that G i7-'^(ri)/M, and u^^v!^~^ in [/ are given. We seek 
for tt'^+^ G such that 

3 
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vl^ ^ \v] \ dx 



v!' W\v}' \ v \ dx 



- [ 2(00 Au^ \v) dx+ [ p''divvdx+ [ 9^+^vids+ [ e^+^usds, (26) 
for ah v£W. Then find q^+^ G H^{Vt)/M. such that 



Vg^+^ I V(/7 ) dx 



di-vvl'^^Lpdx, 



in ^ ' '2St 

for ah if G H\n)/R. The pressure p'^+i G H\n)/R is given by 



p''^'^(pdx= / (g''+^ -Fp^)(/?dx, 
for all if G L^(r2). Finally the velocity u^^^ G {/ is the solution of 

I t;) = J (u^+^ \v^ dx + y (yq''+^ \ dx, 



(27) 



(28) 



(29) 



for ah V G L'^i^f. 



5 



5 Finite element discretization 



A finite element mesh Th of the domain Q into brick elements is considered: Q = Uxerh-f^j where 
each geometric elements K is an hexahedron. For 1 < j < 3, set 

Vh = {^h G C{n); iphlx G Q2, yK G t^] , 
Wh = {v,,Glf; K|n)=Oonr,uri; ^; = Oonro}, 
Uh = [vheV^, {vh\n)={)ond^] , 

The following finite element discretization of the previous projection method is used. 

For /c > 0, assume that E Ph, and u^, u^"^ in Uh are given. Then compute v!^^ E Wh such that 



3 

25t 



n 



I Vh) dx+ [v Vul+^ I Vvh] dx 



2(u;A4|uft) dx+ p'ldiYVhdx+ / (9^+^Wh,i + / el+^Vh,2ds, (30) 
1 ^ ^ Jn JVs Jvs 

for ah E Then find E P/^/M with 

j^(Vql+^\Viph) dx = -^J^d\Yul+^^hdx, (31) 
for all iph £ Ph- The pressure p^^^ E -P/i/M is given by 

p^+Vh,dx= / (g^+i+p^)(^,dx, (32) 
for all ^h ^ Ph- Finally the velocity m^^^ E ?7/t is the solution of 

/ (n^+i I Vh) dx= I (ul+^ I Vh) dx + 5t\ j (Vg^i | Vh) dx, (33) 

for aU f E V^^. 
Remark 



1. The solution of the linear systems (30)-(33) are performed via some preconditioned conjugate 
gradient methods. These methods are well adapted to massively parallel computers. 



2. For equation (31), a coupled preconditionner: diagonal plus optimal conjugate Grahm- 
Schmidt least squares preconditionner (DIAG + LS CGS OPT) is used, see jini dH [T7] . 
Moreover a Lagrange multiplier is used in order to impose the zero mean of Qh^^ , see also [5]. 

3. For the other equations, an incomplete Cholesky ICO preconditionner |14j is sufficient. 
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6 Application to the North Atlantic 



This section is devoted to the numerical simulation of the 3D water circulation in the North Atlantic 
ocean, induced by some mean wind tensions. For this, the methods presented in the previous 
sections are used. A parallel software was developed to solve this kind of problems. The bathymetry, 
and the wind tensions where obtained from the Mercator Ocean Project ( http: / / www.mercator.fr | ) . 
We are grateful to Jean-Marc Molines at INP-Grenoble for his help. 

The model into consideration does not take into account of the temperature, nor the salinity. 
Our aim is, in a first step, to study the currents induced by the wind tensions in the ocean. A 
particular attention is turned on the down, and upwelling near the coast. Our model allows to 
obtain consistent results with reality. It demonstrate that the winds are the main driving forces 
for the global dynamic in the oceans. 

The part of the Atlantic ocean taken into account is delimited in the East by the European and 
the African continents, in the West by the American continent, in the South by the Equator, and 
in the North by the 70° parallel. These two parts lead to some artificial boundary Ti (see figure 
[2|. On this part of the boundary, it is assumed that there is no tangential tensions, and that there 
is no water exchange with outside. 




CD 



Figure 2: Part of the Atlantic Ocean 
The bathymetry is given in figure |3j and|4} 
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Figure 3: Bathymetry of the Atlantic Ocean 




Figure 4: Bathymetry of the Atlantic Ocean, with vertical scale 
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The mean wind during 15 years (1979-1993, ERA15) are represented in figure [5] 
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Figure 6: Mean wind tensions ERA15 (1979-1993) 

The following physical constants are used in the equations: 

• The turbulent cinematic viscosity tensor is = (10*^, 10^, 2.5 • 10^) [mVs], it is based on the 
study of the Reynolds turbulent tensor. 

• The time step is 6t = 2.592 • 10'^ (= 1 month). 

• The final time is set to 100 years 

In the following figures, some numerical results for the currents at different depth are presented. 
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Figure 7: Currents at the surface 




Figure 8: Currents at 500 m depth 
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Figure 13: Currents at 5000 m depth 

The next figures sliow the streamhnes associated to the previous velocity field. The upwelling and 
downwelling are shown near the American and the African coasts. 
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Figure 14: Computed streamlines in the North Atlantic 





Figure 16: Down and upwelling on the American coast 

Let us mention that the travel time for a particle of water, at the surface, to cross the Atlantic 
from the African coast to the South American coast is about 250 days. The travel time of the same 
particle to come back from the South American coast to the African one through the deep water 
is about 12 years. 

Finally a representation of the surface streamlines obtained with our model is represented in figure 



17 This figure can be compared with figure 18 obtained on the web site [I]. 

The main difference between these two figures are in the region of the Caribbean Islands. The 
reasons for these differences are 

• Our model do not take into account the temperature and the salinity. Therefore we do not 
have buoyancy effects. 

• Because of our boundary conditions on the Equator, we could not take into account of the 
South Equatorial current along the East coast of South America. For this we could have to 
compute the currents in the entire Atlantic, but we could not have access to the data for this. 
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Figure 17: Streamlines at the surface 




Figure 18: North Atlantic Gyre [T] 
As already mentioned, our model is able to precisely describe the down, and upwelling near the 
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coast. Moreover it allows to obtain consistent results with measurements, and demonstrate that 
the winds are the main driving forces for the global dynamic in the oceans. 
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